/*==================================================
Project:       Targeting Social Programs
Authors:       Diether W. Beuermann
              Bridget Hoffmann        
              Marco Stampini 
              David L. Vargas
              Diego Vera-Cossio
----------------------------------------------------
Creation Date:    May 2023
Modification Date:   
Do-file version:    01
References:          
Output:             
==================================================*/

/* Regular PMT model using all SISBEN varibles and SISBEN income
This will serve as baseline for comparison		*/


/*==================================================
              0: Program set up
==================================================*/
*Written on STATA 17
drop _all
*set varabbrev off	// no variable abbreviations allowed (personal preference)

/*==================================================
            1: load and transformations
==================================================*/

*----------- 1.0.1 Load CRRA welfare FUN
do "${dir6r}/ados/bstrap_pmt.ado" // load bstap FUN
do "${dir6r}/ados/bstrap_pmt_mh.ado" // load bstap FUN
do "${dir6r}/ados/bstrap_pmt_covfix.ado" // load bstap FUN
do "${dir6r}/ados/baseline_pmt.ado" // load baseline FUN
do "${dir6r}/ados/expanded_pmt.ado" // load expanded FUN

*----------  1.1. Data prep:
use "${dir3r}/01_survey/survey_targeting_r1_long.dta", clear

* new interactions 
cap drop d_loss_illness_inc
cap drop d_reco_ilness_inc
gen d_loss_illness_inc = illness_incidence * d_loss_ilness
gen d_reco_ilness_inc = lag_illness_incidence * d_reco_ilness
lab var d_loss_illness_inc	"Loss: Illness Incidence"
lab var d_reco_ilness_inc	"Recovery: Illness Incidence"


cap drop d_loss_cut_remittance_inc
cap drop d_reco_cut_remittance_inc
gen d_loss_cut_remittance_inc = cut_remittance_incidence * d_loss_cut_remittance
gen d_reco_cut_remittance_inc = lag_cut_remittance_incidence * d_reco_cut_remittance
lab var d_loss_cut_remittance_inc "Loss: Remittances Cut Incidence"
lab var d_reco_cut_remittance_inc "Recovery: Remittances Cut Incidence"

* Globals with shocks
glo loss_rec      d_loss_ilness d_reco_ilness d_loss_cut_remittance d_reco_cut_remittance death divorce bankrupt nat_shock 
glo loss_inc      illness_incidence cut_remittance_incidence death divorce bankrupt nat_shock 
glo loss_rec_inc  d_loss_illness_inc d_reco_ilness_inc d_loss_cut_remittance_inc d_reco_cut_remittance_inc death divorce bankrupt nat_shock 




*Globals with updted survey variables
glo hh_char_srv urban age_feb20_srv prop_kids_srv n_members_srv elementary_srv highschool_srv tertiary_edu_srv children_srv living_couple
glo assets_srv own_washing_machine_srv own_tractor_srv own_moto_srv own_car_srv own_PC_srv own_fridge_srv
glo hh_srv  wall_finished_srv floor_finished_srv cookpower_connected_srv wc_flush_srv kitchen_srv electricity_srv running_water_srv trash_service_srv 
glo labour_srv formal_work_srv informal_work_srv

// Including labour status from SISBEN
glo pmt_1_srv  $hh_char_srv $assets_srv $hh_srv $labour_srv
** make sure all labels are OK 
lab var illness_incidence 			"Illness Incidence"
lab var cut_remittance_incidence 	"Remittances Cut Incidence"
lab var death						"Death"
lab var divorce					"Separation of Spouses"
lab var bankrupt					"Bankruptcy or Business Closure"
lab var nat_shock					"Natural Disaster or Fire"
lab var job_loss                   "Loss: Job lost at any point of the year"
lab var job_loss_formal            "Loss: Job lost formal at any point of the year"
lab var job_loss_informal          "Loss: Job lost informal at any point of the year"

*----------  1.1.2 Estimate PMT for survey data:


//============== Dynamic estimation 

// -------- OLS
local nreps=1000
local base=1.90*1515.1
// 2019 PPP factor https://data.worldbank.org/indicator/PA.NUS.PRVT.PP?locations=CO


baseline_pmt, pmtvars(l_pp_inc_srv $pmt_1_srv) predvars($pmt_1_srv)  reps(`nreps')  base(`base') negonly(no)
mat estmat_static = (1\1\1), r(estt)
scalar basehhtrans = r(estt)[1,10]
scalar coverage = 100*r(estt)[1,2]
baseline_pmt, pmtvars(l_pp_inc_srv $pmt_1_srv)  predvars($labour_srv) negonly($assets_srv) reps(`nreps')  base(`base') fixtrans(basehhtrans)
mat amat = (2\2\2), r(estt)
mat estmat_static = estmat_static \ amat
expanded_pmt,  pmtvars(l_pp_inc_srv $pmt_1_srv) predvars($pmt_1_srv)  reps(`nreps')  base(`base') fixtrans(basehhtrans)
mat amat = (3\3\3), r(estt)
mat estmat_static = estmat_static \ amat
bstrap_pmt,  incvar(d_income_nt) variables(d_gain_work_any_srv d_loss_work_any_srv) reps(`nreps')   base(`base') pmtvars(l_pp_inc_srv $pmt_1_srv) predvars($pmt_1_srv) fixtrans(basehhtrans)
mat amat = (4\4\4), r(estt)
mat estmat_static = estmat_static \ amat

bstrap_pmt_covfix,  incvar(d_income_nt) variables() reps(`nreps')   base(`base') pmtvars(l_pp_inc_srv $pmt_1_srv) predvars($pmt_1_srv) fixtrans(basehhtrans) coverage(10)
mat amat = (5\5\5), r(estt)
mat estmat_static = estmat_static \ amat

bstrap_pmt_covfix,  incvar(d_income_nt) variables(d_gain_work_any_srv d_loss_work_any_srv) reps(`nreps')   base(`base') pmtvars(l_pp_inc_srv $pmt_1_srv) predvars($pmt_1_srv) fixtrans(basehhtrans) coverage(10)
mat amat = (6\6\6), r(estt)
mat estmat_static = estmat_static \ amat 

bstrap_pmt_covfix,  incvar(d_income_nt) variables() reps(`nreps')   base(`base') pmtvars(l_pp_inc_srv $pmt_1_srv) predvars($pmt_1_srv) fixtrans(basehhtrans) coverage(30)
mat amat = (7\7\7), r(estt)
mat estmat_static = estmat_static \ amat

bstrap_pmt_covfix,  incvar(d_income_nt) variables(d_gain_work_any_srv d_loss_work_any_srv) reps(`nreps')   base(`base') pmtvars(l_pp_inc_srv $pmt_1_srv) predvars($pmt_1_srv) fixtrans(basehhtrans) coverage(30)
mat amat = (8\8\8), r(estt)
mat estmat_static = estmat_static \ amat   

bstrap_pmt_covfix,  incvar(d_income_nt) variables() reps(`nreps')   base(`base') pmtvars(l_pp_inc_srv $pmt_1_srv) predvars($pmt_1_srv) fixtrans(basehhtrans) coverage(50)
mat amat = (9\9\9), r(estt)
mat estmat_static = estmat_static \ amat

bstrap_pmt_covfix,  incvar(d_income_nt) variables(d_gain_work_any_srv d_loss_work_any_srv) reps(`nreps')   base(`base') pmtvars(l_pp_inc_srv $pmt_1_srv) predvars($pmt_1_srv) fixtrans(basehhtrans) coverage(50)
mat amat = (10\10\10), r(estt)
mat estmat_static = estmat_static \ amat   

bstrap_pmt_covfix,  incvar(d_income_nt) variables() reps(`nreps')   base(`base') pmtvars(l_pp_inc_srv $pmt_1_srv) predvars($pmt_1_srv) fixtrans(basehhtrans) coverage(70)
mat amat = (11\11\11), r(estt)
mat estmat_static = estmat_static \ amat

bstrap_pmt_covfix,  incvar(d_income_nt) variables(d_gain_work_any_srv d_loss_work_any_srv) reps(`nreps')   base(`base') pmtvars(l_pp_inc_srv $pmt_1_srv) predvars($pmt_1_srv) fixtrans(basehhtrans) coverage(70)
mat amat = (12\12\12), r(estt)
mat estmat_static = estmat_static \ amat    

bstrap_pmt_covfix,  incvar(d_income_nt) variables() reps(`nreps')   base(`base') pmtvars(l_pp_inc_srv $pmt_1_srv) predvars($pmt_1_srv) fixtrans(basehhtrans) coverage(90)
mat amat = (13\13\13), r(estt)
mat estmat_static = estmat_static \ amat

bstrap_pmt_covfix,  incvar(d_income_nt) variables(d_gain_work_any_srv d_loss_work_any_srv) reps(`nreps')   base(`base') pmtvars(l_pp_inc_srv $pmt_1_srv) predvars($pmt_1_srv) fixtrans(basehhtrans) coverage(90)
mat amat = (14\14\14), r(estt)
mat estmat_static = estmat_static \ amat    


/*==================================================
            4: Consolidate data
==================================================*/

** send to stata
cap frame create points
frame change points 

// Static models

drop _all 
svmat estmat_static, names(col)
g model_n = c1 
gen type_m = ""
replace type_m =  "Baseline"  if model_n == 1
replace type_m =  "Updated"   if model_n == 2
replace type_m =  "Expanded"  if model_n == 3
replace type_m = "Dynamic" if model_n == 4
replace type_m = "Baseline-10" if model_n==5
replace type_m = "Dynamic-10" if model_n==6
replace type_m = "Baseline-30" if model_n==7
replace type_m = "Dynamic-30" if model_n==8
replace type_m = "Baseline-50" if model_n==9
replace type_m = "Dynamic-50" if model_n==10
replace type_m = "Baseline-70" if model_n==11
replace type_m = "Dynamic-70" if model_n==12
  replace type_m = "Baseline-90" if model_n==13
replace type_m = "Dynamic-90" if model_n==14
cap drop base _merge 
drop c1



foreach v in covered inclusion_err inclusion_ok exclusion_err exclusion_ok zero epov {
  replace `v' = `v' * 100
replace `v'_uci=`v'_uci*100
replace `v'_lci=`v'_lci*100
}

foreach v in covered hh_benefits inclusion_err inclusion_ok exclusion_err exclusion_ok crra crra_l crra_u zero epov covered_uci covered_lci inclusion_err_uci inclusion_err_lci inclusion_ok_uci inclusion_ok_lci exclusion_err_uci exclusion_err_lci exclusion_ok_uci exclusion_ok_lci crra_uci crra_lci crra_l_uci crra_l_lci crra_u_uci crra_u_lci hh_benefits_uci hh_benefits_lci budget_nat budget_nat_uci budget_nat_lci zero_uci zero_lci epov_uci epov_lci {
  sum `v' if year == 2019 & model_n == 1
  replace `v' = r(mean) if year == 2019 & model_n != 1 & model_n<=4
}

// PPP values 

gen hh_benefits_lcu = hh_benefits
replace hh_benefits = hh_benefits / 1515.1 
replace hh_benefits_uci= hh_benefits_uci/1515.1
replace hh_benefits_lci= hh_benefits_lci/1515.1
// 2019 PPP factor https://data.worldbank.org/indicator/PA.NUS.PRVT.PP?locations=CO

foreach var in budget_nat budget_nat_lci budget_nat_uci {
replace `var'=`var'/ 1515.1 
replace `var'=`var'/1000000
}

save "${dir3r}/03_auxiliar/Budget_estimates.dta" , replace
